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Abstract 

This work presents the formalism and implementation of excited state nuclear forces within 
density functional linear response theory (TDDFT) using a plane wave basis set. An implicit 
differentiation technique is developed for computing nonadiabatic coupling between Kohn-Sham 
molecular orbital wavefunctions as well as gradients of orbital energies which are then used 
to calculate excited state nuclear forces. The algorithm has been implemented in a plane 
wave/pseudopotential code taking into account only a reduced active subspace of molecular 
orbitals. It is demonstrated for the H2 and N2 molecules that the analytical gradients rapidly con- 
verge to the exact forces when the active subspace of molecular orbitals approaches completeness. 
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I. INTRODUCTION 



The past decade has seen time-dependent density functional linear response theory 
(TDDFT)^2aL become the most widely used electronic structure method for calculating ver- 
tical electronic excitation energies 4 ^. Except for certain well-known problem cases such as, 
for instance, charge transferrin and double excitations 9 -, TDDFT excitation energies are 
generally remarkably accurate, typically to within a fraction of an electron Vol t lc>l11112113 . 

Excited state analytical nuclear forces within TDDFT have only been implemented 
recentl y 14 ' 15 ' 16 ' 17 in an attempt to extend the applicability of TDDFT beyond single point 
calculations. One complication has been the fact that TDDFT merely provides excita- 
tion energies, but excited state wavefunctions are not properly defined. The first excited 
state geometry optimization using analytical gradients was presented by van Caillie and 
Amos based on a Handy- Schaefer Z-vector metho d 14 ' 15 . An extended Lagrangian ansatz 
was chosen by Furche and Ahlrichs 16 and Hutted for their Gaussian-type basis set and 
plane wave/pseudopotential implementations, respectively. The latter variant is of particu- 
lar importance for condensed phase applications since it is used in conjunction with periodic 
boundary conditions. In order to ensure completeness, the number of Kohn-Sham (KS) or- 
bitals included in constructing the response matrix in a molecular orbital (MO) basis must 
equal the number of basis functions. Since a plane wave basis typically consists of two 
orders of magnitude more basis functions than a Gaussian-type basis set a complete MO 
formulation of TDDFT is impractical. A solution to this problem is to cast the working 
matix equations directly into a plane wave basis as proposed by Hutter— . Earlier, Doltsinis 
and Sprik 1 ^ have proposed an alternative, active space approach to TDDFT in which only a 
subset of (active) KS orbitals is selected to construct the response matrix. For a large variety 
of excited states, convergence of the corresponding excitation energies has been shown to be 
rapid with respect to the number of orbitals included in the active space^^. In the present 
paper, we shall follow this active space ansatz and derive analytical expressions for excited 
state nuclear forces within an MO basis. In contrast to previous work, we do not rely on a 
Lagrangian formulatio n 16 ' 17 ' 19 , but employ an implicit differentiation scheme instead. This 
has the advantage that we obtain, in addition to the excited state energy gradients, also the 
gradients of KS energies and wavefunctions. The latter may be exploited to compute nona- 
diabatic coupling matrix elements between different electronic states. We have implemented 



the working equations within a plane wave/pseudopotential code^ and we will demonstrate 
numerical accuracy and illustrate the convergence behaviour with respect to the number of 
MOs included in the active space for some prototypical test examples. 



II. THEORY 

The linear response TDDFT eigenvalue problem can be written in Hermitian form as 

Q\F n ) = u 2 n \F n ) , (1) 
where the response matrix f2 is defined as 

Qphcr,p'h'(T' — S a(T i5pp'5hh'(€ p(7 — € ha .) 2 + 2^Je p(T — £hoKpha,p'h'a' \j £pV ~~ € h'a , (2) 

with the coupling matrix 

Kphayh'a* = J dr J dr'^ pa (v)^ ha (r)f^ c (r,v')ijp /a >(r')ij h > a >(v') , (3) 

ip pa and ipho being the KS particle (unoccupied) and hole (occupied) molecular orbitals with 
spin a corresponding to the KS energies e pa and e^, respectively The response kernel 

^'^vh\ + ^-^wBpv> (4) 

containing a Hartree term and an exchange-correlation term is given in the usual adiabatic 
approximation^, i.e. the exchange-correlation contribution is taken to be simply the second 
derivative of the static ground state exchange-correlation energy, i? xc , with respect to the 
spin density p a . 

For the sake of simplicity, the following derivation is for singlet excitations only (extension 
to triplet excitations is straightforward). We shall therefore drop the spin index a. The 
reduced (singlet) response matrix is given by 

Qph,p'h' = Spp'Shh'(^p — £/i) 2 + 4a/ € p — €hKph,p'h'\/ Cp' ~ € h' , (5) 

where 

K ph ,p>h> = J dr y drVp(r)^(r)/ H ,xc(r,r / )V ; p'( r ' / )^'( r ' / ) > ( 6 ) 

and 

/H,xc(r,r') = T ^ r + 5(r-r')|^| . (7) 



Multiplying eq. ((TJ) by (F n | from the left we obtain 



(p„|n|F n > = ul 



(8) 



Differentiation with respect to the nuclear coordinate R a (a = 1, . . . , 3N) for a molecule 
consisting of iV atoms yields 



(9) 



ph p'h' 



where we have used the short-hand notation j^- = f a ; the F^j are the components of 
the linear response eigenvector F n . Carrying out the differentiation of the response matrix, 
eq. © becomes 
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Here we have defined the contracted densities 
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(12) 



r«(r) = ^(r)^-(r) (13) 

In order to compute the excitation energy gradient (eq. (JTUJ)) , we require the nuclear deriva- 
tives of KS orbital energies and wavef unctions, ef and (i — p,h). These can be obtained 
using an implicit differentiation scheme as follows. We start by writing down the KS equa- 
tions in matrix form 



For the full differential of F^ we have 



13 [ dR n 



(14) 



(15) 



jk — JM: 



where = s ^ r y Division by dR a yields 



^ - &v = - E / drH iM r ) = ~ E / dr J dr ' H iA r - r ')C(r') • (16) 

fc k 

On the rhs of eq. we have inserted a delta function, which we now express in terms of 
KS orbitals 

5(r-r') = EV^(rM(r / ) . (17) 
l 

Thus eq. (fTSj) becomes 



where 



= (5ife5y + 8j k 5u)ei + 2n k J ' dr J drT fe /(r)/ H ,xc(r, r^r^r') , (19) 



and 



= y dr^(r)C(r) , (20) 
being the number of electrons occupying orbital k. 

Exploiting the symmetry of the nonadiabatic coupling matrix elements ([20)1. i.e. tpf k = 
—ip^ 1 and therefore ip™ 1 = 0, eq. (JTHj) can be rewritten as 

dm. 



dR a , t 
and for the diagonal terms (i = j) 



E^r ,(*<j) (2i) 



where 

D% = Ef 3 - E\] = (8 u 5 kj + 8 lk Si 3 )(e k - e { ) + 2(n, - n k )K ihlk (23) 
With the definition (|2*3)) eq. (|2~Tf becomes 

= E((V ~~ e h')8p P >5hh' + ^K p , h , )Ph )^l?' (24) 

p'n' 

for particle-hole states, and 
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for all remaining combinations. Eq. (|2"5|) allows us to express the nonadiabatic coupling 
elements between non-particle-hole states analytically as 



+ 4 J2 Ph K aM, 
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The system of linear equations (|24|) is first solved for the particle-hole nonadiabatic coupling 
elements ipp h , which are then inserted into eq. (|26|) to obtain the remaining, non-particle- 
hole, elements. The second term in the numerator of eq. (J26|) is most conveniently evaluated 
by introducing the contracted density 

r 3 (r) = 5> p (r)^(r)C P (27) 

ph 



Then 



J2KijMp h = j dr J C /rT 3 (r)/H, X c(r,r')^(r')^(r / ) = (28) 



W = ^ , (« < J, ij 3 ph) (29) 



Thus eq. (EKJ) becomes 



Similarly the KS orbital energy gradients can now be obtained from the simplified eq. 

£? = m; +4K ■ (30) 

Finally, the nuclear derivative of the KS orbital wavefunction is recovered by unfolding the 
nonadiabatic couplings 

V£(r)=5>(r)C' • (31) 
i 

Equations (1)— (31) have been implemented with periodic boundary conditions using a 
plane wave expansion of the KS MOs at the T point of the Brillouin zone. By making use 
of the periodic boundary conditions, the generalized densities Ti, T 2 , T 3 and can be 
expanded in reciprocal space via the three-dimensional Fourier transform, e.g. 

r 1 (r) = ^r 1 (G)exp(zGr) (32) 

G 

where G is the vector of the recipocal lattice. The Hartree part of the matrix ele- 
ment J dr f dr / ri(r)/H,xc(i", r')r2(r / ) and J dr J dr'Tki{r)f-n tXC (r,r')Tij{r') which enter the 
key equations (fTUj) and (|T9"j) . respectively, can be readily computed in reciprocal space, e.g. 

! dr ! dr' ^(r)— ^-r 2 (r') =n£ ^(G^G) (33) 

II Gf L Q 



whereas the exchange-correlation parts of the matrix elements are calculated via direct 
numerical integration over grid in coordinate space. 



III. TEST RESULTS 

To illustrate the performance and the convergence behaviour of our method, we have 
computed nuclear gradients of the first excited state energies of H 2 and N 2 . The calculations 
were performed using our implementation of the formalism presented here in the CPMD 
package^. All the systems were treated employing periodic boundary conditions and the 
molecular orbitals were expanded in plane waves at the T point of the Brillouin zone. We used 
Troullier-Martins normconserving pseudopotentials^. The excitation energies and nuclear 
gradients were computed within the adiabatic local density approximation^ to the linear 
response exchange-correlation kernel. 

The central idea underlying the present active space approach originates from the ob- 
servation that excitation energies for a large number of electronic transitions exhibit only a 
minor dependence on the size of the response matrix (J2J). This is illustrated in Table 1 for 
the 3cr g — ■> liTg transition of N 2 . A simple two-state HOMO-LUMO response calculation is 
seen to give an excitation energy which is less than 0.2 eV away from an extended treatment 
including all 5 occupied and 100 virtual MOs. Generally, such behaviour is to be expected 
for excitations which can be characterized by only a few low-lying one-electron transitions 
without higher-lying continuum states mixing in. 

In the following, we shall discuss the active space dependence of the excitation energy 
gradients and the KS orbital energy gradients. The upper panel of Fig. 1 displays the 
completeness of the active space as a function of the number of virtual KS orbitals included 
in the space. The integral 



was used as a measure of completeness of the active space. It becomes unity when the active 
space of KS orbitals is complete, i.e. when the total number of the KS orbitals (virtual 
and occupied) equals the number of plane waves used to solve the KS equations. The 
total number of plane waves is 925 for the 6 a.u. cubic box and 40 Ry plane wave cutoff. 
With 450 virtual orbitals included the active space is almost complete and the value of the 




i=l 



(34) 



integral (|H4|) deviates from unity by approximately 10~ 3 which is already comparable with 
the accuracy of the numerical integration. The lower panel of Fig. 1 shows the absolute 
deviation of the analytic derivatives from the respective finite difference values for the the 
first singlet excitation energy, u\, as well as the HOMO and LUMO KS orbital energies, e± 
and €2, of H 2 at a bond length of 1.0 cL.U. clS db function of the number of virtual KS orbitals 
included in the active space. The absolute deviation in analytical gradients vanishes rapidly 
as the number of virtual orbitals is increased and the errors in the analytical gradients of 
different states generally show the same patterns in the dependence upon the number of 
virtual orbitals included in the active space. 

Fig. 2 shows the absolute deviation of the analytical derivative from the respective finite 
difference value as a function of the size of the active space for the first three KS orbital 
energies e, (i = 1, 2, 3) as well as the lowest response matrix eigenvalue LO\ of the N 2 molecule. 
The errors of the analytical gradients are seen to decrease rapidly as the number of orbitals 
included in the active space approaches the number of plane wave basis functions (in this 
case 925 plane waves). For the largest active space the deviations of all energies are of the 
order of 1CT 3 or smaller. At this point, the accuracy of the analytical derivatives is hard 
to assess because the finite difference reference values are also subject to numerical errors. 
We have further checked whether the excited state gradients are invariant under translation. 
The translational contribution to the excitation energy gradient is found to decrease rapidly 
as the active space increses. Interstingly, the underlying ground state calculation exhibits a 
significantly larger translational error than the excitation energy. 

To test the practical value of our derivatives, we have performed geometry optimizations 
of Ni in the first excited state (8 x 5.6 x 5.6 a.u. box, 40 Ry plane-wave cutoff, i.e. 600 basis 
functions). When we include only 100 virtual orbitals in the active space, we obtain a bond 
length of 2.44 a.u., which deviates by 0.02 a.u. from the value of 2.42 a.u. determined by 
a series of single point energy calculations. Upon increasing the number of virtual states to 
200, the optimized bond length comes out as 2.42 a.u., correct to two decimal places. Our 
test calculations illustrate how the size of the active space may be adjusted to achieve any 
desired level of accuracy. For many practical purposes it will be sufficient to work with a 
reduced active space which is significantly smaller than the total number of basis functions. 

The situation is different in the case of molecular dynamics simulations, where the nuclear 
forces need to be essentially exact derivatives of the potential in order to maintain conser- 



vation of energy. Although we have already carried out test excited state MD simulations 
for diatomic molecules, those results do not offer any additional insight into the general 
performance and convergence pattern of our method. We have not yet applied the formal- 
ism presented here to perform more realistic excited state MD simulations of polyatomic 
molecules , because our current implementation does not yet make use of more efficient 
iterative techniques, such as Lanczos algorithm or related schemes^, to solve the response 
eigenvalue problem These numerical techniques, however, are standard and we plan to 
exploit them in future implementations. The scope of this article is merely the presentation 
of the formalism and the analysis of the convergence behaviour with respect to the choice 
of the active space. 

We would like to emphasize, however, that the method described here is capable of 
providing additional information beyond excited state energy gradients. Fig. 3 shows, for 
instance, the nonadiabatic coupling strength between the second and third KS orbitals, -02 
and ^3, for the H 2 molecule as a function of its bond length. The nonadiabatic coupling 
values obtained from eqn (J29j) exhibit a singularity at the crossing point between the two KS 
orbital energies, as one would expect due to the KS energy difference in the denominator. 
This feature of our formalism may be exploited in future applications of TDDFT beyond 
the Born-Oppenheimer approximation. 

IV. CONCLUSIONS 

We have developed and implemented an novel, alternative formalism to calculate analyti- 
cal nuclear forces for TDDFT excited states within a plane wave/pseudopotential framework. 
In addition to the excited state energy gradient, our method also provides molecular orbital 
wavefunction as well as energy derivatives at a small computational overhead compared to 
the vertical excitation energies. The latter may, for instance be employed as a powerful tool 
for understanding and interpreting various chemical phenomena such as molecular structures 
and reactivities^. A fundamental quantity in the present formalism are the nonadiabatic 
coupling elements in the molecular orbital basis which are obtained as direct solutions of a 
system of linear equations. These matrix elements may provide the basis for the calculation 
of nonadiabatic couplings between the many-electron adiabatic wavef unctions. Trial calcu- 
lations on the prototypical test molecules H 2 and N 2 demonstrate that our implementation 



reproduces the exact gradients when the number of molecular orbital included in the active 
space approaches the number of plane wave basis functions. Excited state geometry opti- 
mization of N 2 using different active spaces show that for many practical purposes it will be 
sufficient to work within a relatively small active space. The size of the latter may be tuned 
to achieve any desired level of accuracy. 
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APPENDIX A: CALCULATION OF ^ IN PLANE WAVES 

The matrix elements of the derivatives of the KS Hamiltonian are required for the cal- 
culations of the KS molecular orbitals gradients. Since the kinetic and exchange-correlation 
energies do not depend directly upon the atomic positions, the matrix element of the deriva- 
tive becomes (local/nonlocal pseudopotential and electrostatic interaction contributions): 

U±±ij ^ U jjppjocal , u tjPP , nonlocal . u ZJ es ( A 1 \ 

All these matrix elements are computed in reciprocal space. The matrix element of the 
derivative of the local pseudopotential has the following form 

JLh»*»* = -n y £iGV£ ad {G)S I {G)T iJ (G) (A2) 

where Rj denotes the atomic position and the structure factor Sj = exp(zGRi) of nuclei 
I, Tij(G) is the three-dimensional Fourier transform of the contracted density (f 13)1 . The 
nonlocal pseudopotential contribution is 







tjPP , nonlocal \ ^ 



fj E 1 " Pi t? v 



(A3) 



where the contribution from the projector derivative -t-^j is calculated in the standard way2£. 
The contribution from the electrostatic energy is computed in the form 

' ; ff|-^4 r l'( G K(G)5i(G) (A4) 



1 G^O 

where n^iG) is the gaussian core charge distribution for nuclei / in reciprocal space. 



APPENDIX B: THIRD FUNCTIONAL DERIVATIVES OF THE EXCHANGE- 
CORRELATION FUNCTIONAL 

1. Third derivative of Vosko-Wilk-Nusair correlation 
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with x = y/r s , r s - y ^ pi 
A = 0.0310907, b = 3.72744, c = 12.9352, x = -0.10498. 
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Using relation (jB4|) . eqn ()B7|) can be rewritten as 
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and thus eqn (jB6|) becomes 
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One arrives at eqn (|B14|) by substituting eqn (|B9j) into eqn ()B13j) . For instance, 
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2. Third derivative of Slater exchange 



E x = j pe x dr , (B19) 

where 

e x = Caph , C = - jj 0V , a = (B20) 
The third derivatives can be readily computed 

= -Ca^p-"- 1 (B21) 
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active space 


AE KS lo/lv 5o/lv 5o/5v 5o/50v 5o/100v 


exc. energy 


8.39 9.47 9.40 9.40 9.33 9.29 



TABLE I: Dependence of N2 TDLDA excitation energy ^Hg, 3a g — ► lir g ) in eV using plane waves 
(p.w.) with a 70 Ry cutoff in a 10 ao periodic box on the number of occupied (o) and virtual (v) 
Kohn-Sham orbitals included in the active space. AEks is the unperturbed Kohn-Sham energy 
difference. 
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FIG. 1: (a) The integral (|34j) as a function of the number of virtual orbitals included in the active 
space, (b) Absolute deviation of finite difference and analytic derivatives of the Kohn-Sham HOMO 
and LUMO energies of H2 at a bond length of 1.0 a.u. as a function of the number of virtual Kohn- 
Sham orbitals included in the active space, (c) Absolute deviation of finite difference and analytic 
derivatives of the first singlet excitation energy of H2 at a bond length of 1.0 a.u. as a function 
of the number of virtual Kohn-Sham orbitals included in the active space. The calculations were 
carried out in a cubix box of length 6 a.u. with periodic boundary conditions and a plane wave 
cutoff of 40 Ry. 




FIG. 2: Absolute deviation of the analytic derivatives from the respective finite difference values 
for the the first singlet excitation energy, lo±, as well as the three lowest KS orbital energies, e« 
(i = 1, 2, 3), of N2 at a bond length of 2.0 a.u. as a function of the number of virtual Kohn-Sham 
orbitals included in the active space. The calculations were carried out in a cubix box of length 6 
a.u. with periodic boundary conditions and a plane wave cutoff of 40 Ry. 
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FIG. 3: Upper panel: Absolute value of the nonadiabatic coupling matrix element between the KS 
orbitals ip2 and ip3 along the molecular axis of H2 as a function of the bond length. Lower panel: 
KS orbital energies, €2 and £3 of H2 as a function of bond length. The calculations were carried 
out in a periodic orthorhombic box of size 8 x 5.6 x 5.6 a.u. 3 using a plane wave cutoff of 40 Ry. 



